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Spectroscopic modes provide the most sensitive probe of the very weak interactions responsible for 
the properties of the long- wavelength cycloid in the multiferroic phase of BiFeOs below Tn ~ 640 K. 
Three of the four modes measured by THz and Raman spectroscopies were recently identified using 
a simple microscopic model. While a Dzyaloshinskii-Moriya (DM) interaction D along [—1,2,-1] 
induces the cycloid with wavevector (27r/a)(0.5 + ^,0.5,0.5 — 5) {5 ^ 0.0045), easy-axis anisotropy 
K along the [1,1,1] direction of the electric polarization P induces higher harmonics of the cycloid, 
which split the modes at 2.49 and 2.67 meV and activate the $2 mode at 3.38 meV. How- 
ever, that model could not explain the observed low- frequency mode at about 2.17 meV. We now 
demonstrate that an additional DM interaction D' along [1,1,1] not only produces the observed 
weak ferromagnetic moment of the high- field phase above 18 T but also activates the spectroscopic 
matrix elements of the nearly- degenerate, low-frequency ^0 and $1 modes, although their scatter- 
ing intensities remain extremely weak. Even in the absence of easy- axis anisotropy, D' produces 
cycloidal harmonics that split ^1 and activate $2- However, the observed mode frequencies and se- 
lection rules require that both D' and K are nonzero. This work also resolves an earlier disagreement 
between spectroscopic and inelastic neutron-scattering measurements. 

PACS numbers: 75.25.-j, 75.30.Ds, 78.30.-j, 75.50.Ee 



I. INTRODUCTION 

As the only known room- temperature multiferroic, 
BiFeOs continues to attract a great deal of attention. 
Multiferroic materials offer the tantalizing prospect of 
controlling magnetic properties with electric fields or 
electric polarizations with magnetic fields^. Although 
the ferroelectric transition temperature^ Tc ^ 1100 K 
of BiFeOs is far higher than its Neel temperature^"— 
Tn ^ 640 K, the electric polarization P is enhanced by 
its coupling to the long- wavelength cycloid below Tn \&\ . 
As a result, the magnetic domain distribution below Tn 
can be manipulated by an electric fi.e\3^^^. 

Before BiFeOs can be used in technological applica- 
tions, however, it is essential to understand the micro- 
scopic mechanisms and interactions responsible for its 
magnetic behavior. At frequencies above a few meV up to 
about 70 meV, the spin-wave (SW) spectrum of BiFeOs 
has been used^i^ to determine the nearest-neighbor and 
next-nearest neighbor exchange interactions Ji ~ —4.5 
meV and J2 ^ -0.2 meV between the S = 5/2 Fe^+ 
spinsi^ on a pseudo-cubic lattice with lattice constant 
a ~ 3.96 A. As shown in Fig. 1(a), Ji is the antiferro- 
magnetic (AF) interaction between spins on neighboring 
(1,1,1) planes separated by c = a/>/3 while J2 is the AF 
interaction between neighboring spins on each hexagonal 
layer. 

Below Tn, a long- wavelength cycloid with wavevector 
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;= [-1,2,-1] 



Q = (27r/a) [0.5 + (5,0.5,0.5 - (^] {S ^ 0.0045 
produced by the Dzyaloshinskii-Moriya (DM) interaction 
D = Dy' along = [—1,2,-1] (all unit vectors are 



x' = [1,0,-1] 



FIG. 1: (Color online) (a) The pseudo-cubic cell with ex- 
change interactions Ji and J2 as well as the polarization direc- 
tion z cutting through two hexagonal planes, (b) For domain 
1, a schematic of the spins along the x' axis showing their ro- 
tation about y^ Due to the DM interaction = D'z\ spins 
rotate by r about z in the x^y^ plane. 



assumed normalized to one). As shown in Fig. 1(b), the 
spins of the cycloid lie predominantly in the (-1,2,-1) 
plane normal to y^ 

Whereas the high-frequency portion of the SW spec- 
trum determines the Heisenberg exchange interactions, 
the low- frequency modes measured by THz^^^^^ and 
Ramani^"— spectroscopies can be used to determine the 
small microscopic interactions that control the cycloid. 
Four modes have been detected at frequencies^^ of 2.17, 
2.49, 2.67, and 3.35 meV. By comparison, a model with 
the single DM interaction D only produced a single 
spectroscopically-active mode labeled ^1 at about 2.37 
meV. 

A more realistic mode l^^i^^ also contains the easy-axis 
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anisotropy K along = [1, 1, 1], parallel to the electric 
polarization P. When > 0, splits into two and $2 
at 3.38 meV is activatedi^. Although this model success- 
fully described the upper three spectroscopic modes, with 
predicted frequencies very close to the measured frequen- 
cies, it failed to explain the low- frequency 2.17 mode. In 
addition, it provides conflicting estimates for K based on 
spectroscopic and inelastic neutron-scattering measure- 
ments. 

Several author^— have examined the effects of an- 
other DM interaction = D'z' between neighboring 
hexagonal layers. For a G-type AF, D' produces a weak 
ferromagnetic moment along due to the canting of the 
uniform moments on each hexagonal plane. The moment 
Mo = 2/iB'S'oy' ^ 0.03/iBy' was subsequently observed in 
the metamagnetic phase^i^^ above 18 T. Below 18 T, D' 
was predicted^^ to induce an oscillatory component of 
the cycloid along y\ which has recently been confirmed 
by neutron-scattering measurements^^. 

Based on a model that includes both D and D' in addi- 
tion to the easy-axis anisotropy we evaluate the spin 
state and spectroscopic modes of BiFeOa. Even when 
K = 0^ D' induces higher harmonics of the cycloid that 
split ^1 and activate $2- More remarkably, D' activates 
^0 and <l>i at the cycloidal wave vector. 

We believe that these nearly-degenerate modes are re- 
sponsible for the low- frequency 2.17 meV peak observed 
in spectroscopy measurements. Although a model with 
K = {) can produce four spectroscopic modes, the ^1 se- 
lection rules are reversed and their mode frequencies are 
too small. Therefore, both D' and K are required to ex- 
plain the experimental measurements. With D' ^ 0.054 
meV, corresponding to the observed value^^^^ So = 0.015, 
we estimate that D ^ 0.11 meV and K ^ 0.0035 
meV, which also provide a good description of inelastic 
neutron-scattering measurements^ below 5 meV. 

This paper is divided into seven sections. Section II 
constructs the spin state of BiFeOa. Section III evalu- 
ates the spin dynamics of that state. Section IV evalu- 
ates the spectroscopic modes of that state, and Section 

V discusses the selection rules for those modes. Section 

VI discusses the inelastic neutron-scattering spectrum for 
the low- frequency modes. Section VII contains a brief 
conclusion. Results for the SW intensities are provided 
in Appendix A. The polarization and magnetic matrix 
elements are provided in Appendix B. 



II. SPIN STATE 

With P = Pz\ the three magnetic domains have 
cycloidal wavevectors Q = (27r/a)[0.5 + (5,0.5,0.5 — S] 
(domain 1), (27r/a)[0.5, 0.5 + (5,0.5 - 6] (domain 2), or 
(27r/a)[0.5 + (5, 0.5 — (5, 0.5] (domain 3). By contrast, the 
G-type AF stabilized by a magnetic field^i^^, doping^!, 
or in thin films^^ has wavevector (27r/a)[0.5, 0.5, 0.5]. In 
our discussion of the selection rules governing the spec- 
troscopic modes in Section V, we will assume that all 



three domains are equally populated. Since the spin state 
and dynamics are the same for all three domains, we 
now concentrate on domain 1 with = [1,0,-1] and 

= [-1, 2, —1], as shown in Fig. 1(b). 

The spin state and SW excitations of BiFeOs are eval- 
uated from the Hamiltonian 

H = — Ji ^ Si • Sj — J2 ^ Si • Sj — K^^Siz''^ 



id) 



^ — ^Rj=Ri+a(x— z) 

-^'Er (-1)^"'/^z'.(5,xS,). 

^ — 'Rj=Ri+ax,ay,az 

The first and second exchange terms contain sums (i, j) 
and (i, jy over nearest and next-nearest neighbors on the 
pseudo-cubic lattice. The third term arises from the easy- 
axis anisotropy along z' and the fourth term from the DM 
interaction with D = Dy' . 

Compared to the model for BiFeOs introduced in 
Ref.[2^ and studied in our earlier work^^, H adds the 
DM interaction = D'z' . This term alternates in sign 
with increasing z': (— 1)^^^'/^ changes sign from layer n 
to layer n+1 so the DM interaction (^ — l^^iz'/^JJ' between 
layers n and n+1 has opposite sign to the DM interaction 
between layers n + 1 and n + 2. Hence, the DM interac- 
tion has the same wavevector (27r/a)[0.5, 0.5, 0.5] as a 
G-type AF. 

Because S ~ 1/222, a unit cell containing M = 222 
sites within each of two neighboring (1,1,1) planes is used 
to characterize the distorted cycloid. In zero magnetic 
field, the cycloid can be expanded in odd harmonics^i^ 
of the fundamental wavevector Q (even harmonics are 
also required in non-zero fields). If Sy' (R) is proportional 
to 6'a;/(R), then 



S,>{K) = (-1)«.'/V1-kV52-^.'(R)' 

sgn(sin(27ri5i?x'/<^))' (2) 
Sy'{R) = K^S^ - 5^KR-)^sgn(sin(27r(5i?xVa)), (3) 

CO 

C2.n+icos((2m+ l)2^(5i?^Va). (4) 

m=0 

Odd-order coefficients C2m+i in Sz'CR) satisfy 
Z)m=o^2m+i = 1. Although Sy'^R) (unlike 6'a,/(R) and 
^'^'(R)) does not change sign from one layer to the next, 
the average value of Syf{Il) vanishes and there is no 
net moment in any direction. The ratio Syf{Il)/Sx'(R) 
has magnitude /^/Vl — which is proportional to 
\D' / Ji\ <C 1. Hence, the tilting angle r indicated in 
Fig. 1(b) satisfies the relation tanr = kI\/\ — ^ 
Although the cycloid remains coplanar for each hexag- 
onal layer, the cycloidal planes rotate by 2r from one 
layer to the next. 

The parameters of the spin state are evaluated by min- 
imizing the energy E = {H) in a unit cell x'y' z' of di- 
mensions 15,000a X a X 2c containing two (1,1,1) layers. 
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Open boundary conditions are employed along the x' di- 
rection. With the exchange interactions J\ = —4.5 meV 
and J2 = —0.2 meV fixed at the values required to de- 
scribe the SW spectrum^'^ at high frequencies, the four 
variational parameters are (5, a^, C3, and C5. A solution 
with S = 1/222 is obtained by varying the DM interac- 
tion D for fixed K. After minimizing the energy, we ver- 
ify that the corresponding spin state provides at least a 
metastable minimum by checking that the classical forces 
on each spin vanish. 

With a magnetic field oriented along z', the metamag- 
netic state observed^i^^ above 18 T can be written 

51 = 5 (cos ^ cos (/), cos ^ sin (/), sin ^) , (5) 

52 = 5'(— cos6>cos0, cos^sin(/), sin^), (6) 

for Rz' = 2mc and (2m + l)c, respectively. Extrapolat- 
ing to zero field with ^ = 0, we obtain tan 2^ = D'/Ji. 
Hence, the weak ferromagnetic moment of the metamag- 
netic phase is 



0,11 



Mo = 2/ibSo = 2/iB^sin( 



Ji ' 



(7) 



independent of and J 2- Using Ji = —4.5 meV 

and the experimental resultj^i^^ 5*0 = 0.015, we estimate 
that \D'\ = 0.054 meV, which is is slightly larger than 
the estimate \D'\= 0.046 meV provided in Ref.[24]. 

For the distorted cycloid given by Eqs.(|2]|4]), it is 
straightforward to show that if (5 <C 1, then n ^ D' /2Ji. 
Therefore, the maximum cycloidal spin |5'2y/(R)| equals 
the weak ferromagnetic spin Sq of the metamagnetic 
phase. For the tilting angle, we estimate r ~ 0.34°, a 
bit smaller than the recent neutron-scattering^^ estimate 
of - 1°. 

In Fig. 2 (a), we plot the DM interaction D versus Sq 
for several values of the anisotropy K ranging from to 
0.0035 meV. For if = and 0.0005 meV, D increases 
slightly with 5*0. But for K > 0.001 meV, D decreases 
with Sq. Nevertheless, the variation of D with So is 
rather modest. 

By contrast, the higher harmonics of the cycloid ex- 
hibit a much stronger variation with Sq. Fig. 2(b) reveals 
that the ratio C3/C1 increases with for all K. Since 
Ci = 1 - E„=i C2„+i and ICsl « IC3I, Ci « 1 - C3 and 
C3/C1 w Call + C3). For X = and So > 0, C3 > and 



n=0 



2„+i)'«i(l-2C3) <^. (8) 



Because the interaction energy is optimized when the 
spins lie in the x'y' plane, higher harmonics favor the 
z' nodal regions of the cycloid. When Sq is sufficiently 
smah and > 0, C3 < and (Sf^,) > 1/2 so that higher 
harmonics favor the antinodal regions of the cycloid. 
Experimentally, the ratio of the neutron-scattering in- 
tensity from the third to the first harmonics is given by 

Notice that the third (and higher) harmonics can van- 
ish for nonzero 5*0 and K. When 5'o = 0.015, C3 < 



0.107 



r 

0.0035 meV 


r 

(a) 


0.003 meV 




0.0025 meV ^ 




_0.002 meV 




'0.0015 meV 
0.001 meV 

n nnnc; ma\i ' 




K=:0 meV 

1 1 



0.03 

0.02 
0.01 


_ -0.01 
-0.02 



o 



-0.03 
-0.04 
-0.05 
-0.06 
-0.07 



1 




K = meV___-- — — ^ 




0.0005 meV^ ■ — ^ 




0.001 meV__ — ■ — ■ — ^ 




j.0015 meY,^ ' ^ 




0.002 meV • — — 

0.0025 meV_^^ • ^ 


^ " 


0.003 meV 




"0.0035 meV • — ' 

1 





0.005 



0.01 



0.015 



FIG. 2: (Color online) (a) The DM interaction D and (b) the 
ratio of harmonics C3/C1 versus So for several values of K. 



when K is less than about 0.001 meV and C3 > when K 
is greater than about 0.001 meV. For i^T ^ 0.001 meV, the 
higher harmonics of the cycloid vanish and {Siz''^) = 1/2. 



III. SW EXCITATIONS 

The SW frequencies are calculated using the equations- 
of-motion technique for non-collinear spins outlined in 
Ref.[3l|. A unit cell containing M = 222 sites on each of 
two hexagonal layers is constructed to evaluate the 2M 
SW frequencies a;n(q). SW intensities are obtained from 
the spin-spin correlation function defined by Eq. (|A9p in 
Appendix A. In the absence of damping, the inelastic 
scattering cross section 6'(q, a;) can be expanded as the 
sum over delta functions at each frequency: 

5(q,u;) = ^(l-(fe/q)2)<5(u;-u;„(q))ej(q). (9) 

n,Q! 

The amplitudes S^a{ci) are evaluated using Eq. ()Alip . 



4 



K = 0.001 meV 



K = 0.002 meV 




3.5 



FIG. 3: (Color online) The SW modes of BiFeOs versus r]/5 
for wavevector (27r/a)(0.5 + ?7, 0.5, 0.5 — 77). Dashed lines 
show all possible excitations and the solid lines show only 
those modes with significant intensity above a threshold value. 
All three plots take Sq = 0.015 and D' = 0.054 meV. $0 
(black dot) has a very large y' MR matrix element. The 
low-frequency mode (brown dot) has both ^0 and ^^^^ con- 
tributions with nonzero x' and y' MR matrix elements, re- 
spectively. Whereas ^^2^ (red) has nonzero y' MR matrix 
element, ^f'^^^ (blue) and ^f'^^^ (green) have nonzero x' and z' 
matrix elements, respectively. The EM mode with component 
y' coincides with 



For fixed S'o = 0.015, the SW frequencies are plotted 
in Fig.3 for K = 0.001, and 0.002 meV. Although 
there are 2M modes for every wavevector 27r/a(0.5 + 
77, 0.5, 0.5 — 77), plotted by the dashed lines, only a few of 
those modes have any significant intensity. Modes with 
intensity above an arbitrary cutoff are plotted in the dark 
lines. 

When K ^ 0.001 meV in Fig.3(b), the higher har- 
monics of the cycloid vanish and the SW frequencies 
are similar to those for S'o = and K = discussed 
in Ref.[l9[. In the absence of harmonics, de Sousa and 
Moored labeled the SW frequencies uOnijnQ) (n = 1 or 
2) of a one-dimensional cycloid at multiples m of the cy- 
cloidal wavevector Q = 27r^/a as <l>m and ^rn- Using 
an extended zone scheme and assuming that \m\5 <C 1, 
uJn{mQ) can be approximated by <l>^ = <l>i|m| and 
= ^iVT+TT?. Thcsc relations imply that <l>i = ^0, 
as seen in Fig. 3(b), and that the ^±m and ^±m modes 
cross without repulsion at the zone center q = Q and 
zone boundary q = 0. 

Whether produced by the tilt r or by the anisotropy 
higher odd harmonics of the cycloid introduce higher 
even harmonics in the Hamiltonian H. A 2mQ poten- 
tial will split the ^±m and ^±m modes. As shown in 
Figs. 3(a) and (c), the new m = 1 eigenmodes are labeled 



and 



Notice that ^0 and ^^"^ are nearly 
degenerate for all K. Although too small to see in Fig.3, 
even ^±2 are split by anharmonicity. 



IV. SPECTROSCOPIC MODES 

Because the wavelength of far infrared light greatly 
exceeds atomic length scales, the SW modes measured 
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FIG. 4: (Color online) The evolution of the predicted modes 
with anisotropy K taking So = 0.015 and D' = 0.054 
meV. The horizontal dashed lines are the spectroscopic mode 
frequencies^^. 



by THz and Raman spectroscopies lie at the zone center 
q = Q or 77 = (5. A magnetic resonance (MR) mode has 
nonzero matrix element ((5|Mq,|0), where |0) is the ground 
state and |^) is an excited state with a single magnon 
of wavevector Q. An electromagnon (EM) mode has 
nonzero matrix element {6\P^^\0) so that the induced 
polarization directly couples the ground state to the ex- 
cited state. 

In order to evaluate the MR and EM matrix elements, 
we must first express the magnetic moment M and in- 
duced polarization P^"^ operators in terms of the spin op- 
erators Si. The magnetic moment M = 2/iB XIr con- 
tains a sum over the 2M unique sublattices. In BiFeOs, 
the coupling between the cycloid and electric polariza- 
tion is produced by the inverse DM mechanisnt^ — with 
induced polarization 



>ind 



where the sum is restricted to the 2M sublattices us- 
ing periodic boundary conditions. Within each (1,1,1) 
plane, e^j = ^/2ax.' connects spins at sites and Hj. 
So if {0\SiX Sj\0) points along y', then (0|P^"^|0) points 
along z^ 

Expressions for the matrix elements {S\Ma\0) and 
{S\P^^\0) are provided in Appendix B. Although there 
is no simple relation between the MR matrix elements 
and the SW intensities, the MR and EM modes only 
appear at mode frequencies n with S'^l^,{S) > 0. Gen- 
erally, <l>n modes with {6\My'\0) 7^ also have nonzero 

SW intensities S^l,{6) and sl?l,{S). Hence, those modes 
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excite spins within the x' z' plane of the cycloid (neglect- 
ing its small tilt). On the other hand, modes with 

((5|M^/|0) 7^ or ((5|M^/|0) 7^ also have > 0. 

Hence, those modes excite spins out of the x' z' plane. 

Zone-center modes with nonzero MR matrix elements 
are indicated by the filled circles in Fig. 3. In addition to 
having an enormous SW intensity, the "zero" -frequency^ 
^0 mode has a very large MR matrix element (for K = 
0.0035 meV and = 0.015, \{6\My^\0) \ ^ 8400/iB). The 
2Q potential splits the degenerate modes into ^^^^ 
(((5|M^/|0) ^ 0) and {{S\M,^\0) + 0). The EM 

(((^|P^^^|0) 7^ 0) always coincides with ^^^^ Similarly, 
the smaller 4Q potential splits the $±2 modes. Due to 
its hybridization with <l>o, ^^2^ becomes spectroscopically 
active with ((5|M^/|0) 7^ 0. 

The predicted mode frequencies are plotted versus 
anisotropy for 6*0 = 0.015 in Fig. 4. Both ^\ ' ' and 



(1,2) 



cross near K = 0.001 meV. At 7^ = (5, <l>i has 
no SW intensity and is not spectroscopically active. But 
at 7^ = 0, this mode is responsible for important features 
in the inelastic-scattering spectrum discussed in Section 
VI. 

For K = 0.0035 meV, the mode frequencies are plotted 

versus in Fig. 5 (a), where D and D' are evaluated in 

terms of So for fixed S = 1/222. While the predicted 

spectroscopic mode frequencies decrease slightly with 5*0, 
(2) 

^\ slightly increases. 

When 6*0 = 0, the ^0 and <l>^^^ modes at the zone cen- 
ter T] = S have no SW intensity and their MR matrix el- 
ements vanish. But when /So > 0, the DM interaction 
with wavevector (27r/a)[0.5, 0.5, 0.5] hybridizes ^0 with 
^^^'^'^ and with <l>o. Consequently, their MR matrix 
elements become significant. 

In Fig. 5(b), the mode frequencies and MR matrix el- 
ements of ^0 and ^'i^ are plotted versus together 
with the very small SW intensities of those modes for 
K = 0.0035 meV. As expected from perturbation the- 
ory, the matrix elements {S\Ma\0) grow linearly with 
^0 ~ \D' / Ji\. Moreover, they scale like the square root 
of the SW intensities Sa'a'i^)- Therefore, these modes 
are both spectroscopically and dynamically activated by 
the tilt of the cycloid. It is remarkable that the MR ma- 
trix elements of ^0 and become so large while their 
SW intensities remain extremely weak. 

The dashed horizontal lines in Fig. 4 correspond to the 
four measured spectroscopic frequencies of BiFeOs. We 
believe that the nearly-degenerate ^0 and <I>^^^ modes are 
responsible for the observed low- frequency peak at 2.17 
meV. Recall that those two modes only appear when the 
cycloid is tilted away from the x'z^ plane by the DM 
interaction along z^ The best overall fit to the ob- 
served mode spectrum is obtained with K ^ 0.0035 meV. 
MeasurediS and predicted mode frequencies are summa- 
rized in Table I. 

With So = 0.015 and K = 0.0035 meV, the harmon- 
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FIG. 5: (Color online) (a) The frequencies of the predicted 
modes versus 6*0 for K = 0.0035 meV. Horizontal dashed 
lines are the measured spectroscopic frequencies^, (b) The 
MR matrix elements \{d\Mc,\0)\/ /j.b for ^0 (solid) and <I>^^^ 
(dashed) versus 6*0 for K — 0.0035 meV. Also plotted are the 
intensities Sa'a' (S) of those modes (a — y' for a — x' 
and a = x or z for ol — y') with M = 222. The normalized 
matrix element \{6\Ma\0)\/ fiB 3^'^' i^V^'^ is independent of 
So. The dash-dot curve plots the MR matrix element for 
with a — y' . 



ics of the cycloid have the ratio C3/C1 = —0.050 or 
(Ci/Ca)^ = 400. Elastic neutron-scatteringii and NMR 
measurements^^ indicate that {Ci/C^Y ^5, 
respectively. However, the NMR measurement may over- 
estimate the third harmonic due to the high ^^Fe isotope 
content of the sample^^^. Our estimate for (Cs/Ci)^ is in 
very good agreement with the elastic neutron-scattering 
result. 
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V. SELECTION RULES 

We now consider the selection rules for the THz 
modesi^ii^ for a sample with the single polarization do- 
main P = Pz', where = [1,1,1]. As mentioned in 
Section II, the three possible magnetic domains have 
wavevectors (27r/a)(0.5 + 0.5, 0.5 - (5), (27r/a)(0.5, 0.5 + 
(5, 0.5 - S), and (27r/a)(0.5 + 0.5 - S, 0.5). Since these 
domains have the same energy, we expect them to be 
equally populated. The mode spectrum was measured 
for crossed fields hi = [1,-1,0] and h.2 = [1,1,0]. 

To predict the selection rules for BiFeOs, hi and h.2 
are expressed in terms of the cycloidal unit vectors x', 
y', and W as 

hi = (x' - V^yO/2, 

h2 = x72 + V3y/6 + V2/3z', (11) 
in domain 1 with x' = [1, 0, —1] and y' = [—1, 2, —1]; 

hi = -(x' + \/3y')/2, 

h2 = x'/2 - V^y'/Q + a/273z', (12) 

in domain 2 with x' = [0, 1, —1] and y' = [—2, 1, 1]; and 

hi = x', 

h2 = (y' + \/2z')/\/3, (13) 

in domain 3 with x' = [1, —1,0] and y' = [1, 1, —2]. Al- 
though the following discussion assumes that all three do- 
mains are equally populated, our qualitative conclusions 
remain unchanged even if one or two domain populations 
dominate the sample. 

While ^[^^ ((^|M^/|0) ^ 0) and <I>^^^ ((^^|M^/|0) ^ 0) 

should appear in both fields hi and h2, ^i {{S\Mz' |0) 7^ 
0) should only appear in field h2, which contains a z' 
component. This agrees with the selection rule observed 
by Talbayev et al. [l^. But Nagel at al. [l5[ recently 

found that ^\ ^ survives in field h2, although with dras- 
tically reduced intensity. Notice that the position of ^1^^ 

above ^\ ^ requires that K > 0.001 meV. Therefore, both 
nonzero K and are required to explain the spectro- 
scopic frequencies and selection rules. 

Whereas Talbayev et ali^ found that the low- 
frequency mode appears only in field hi, our model indi- 
cates that the nearly-degenerate ^0 {{^\Mx'\0) ^ 0) and 
^[^^ {{S\My^\0) 7^ 0) modes should appear in both fields 
hi and h2. However, more precise THz measurements^^ 
have recently detected the low-frequency mode in both 
fields hi and h2. At 4 K, Nagel et ali^ even observed 
distinct low- frequency peaks at 2.03 and 2.26 meV. The 
observed three- fold splitting of the 2.03 meV peak in a 
magnetic field may help to distinguish ^0 and 

To address the observability of the THz modes more 
carefully, we evaluate the spectroscopic intensities /(hi) 
and /(h2) for each mode. The spectroscopic intensity for 



TABLE L Spectroscopic Frequencies, Matrix Ele- 
ments, and Intensities 















Measured u (meV) 


2.17 


2.49 


2.67 


3.38 


Predicted cj (meV) 


2.03/2.05 


2.53 


2.75 


3.40 


MR index a 










y' 


K5|M„|0)|/mb 




2.50/1.86 


3.96 


4.59 


1.01 


mPy,\o)\/\ 










12.2 





Intensity index 


a' 








x\z' 






4.94 X 10"V 


19.7 


18.1 


5.43, 






3.05 X 10"^ 






2.35 


/(hi)/Ml 
/(h2)/Ml 




4.75 





10.54 


0.51 




1.58 


10.47 


3.51 


0.17 



any mode is given b}^ 

i{h) = j2himM^\o)\\ (14) 

a 

Averaging over the three domains, we find 

I{h,) = \{\{5\M^,\0)\' + \{5\My,\0)\^], (15) 
I{h2) = ^{\{5\M^,\0)\' + \{5\My,\0)\'} 

+ lmM^'\0)\\ (16) 

For {S\Ma\0) ^ 0, /(hi)//(h2) = 3 for any mode (like 
and $J^^) with a = a;' or y' while 

(2) 

/(hi)//(h2) = for any mode (like ^i ) with a = z' . 

The spectroscopic intensities for K = 0.0035 meV and 
^0 = 0.015 are summarized in Table I. These numer- 
ical results indicate that ^i^^ and ^i^^ should be the 
strongest of the four modes, in agreement with the THz 
results^^'^^. Surprisingly, Table I indicates that the in- 
tensity /(h2) of <I>2^^ is roughly 20 times smaller than 
that of By contrast, recent THz measurement^ 

indicate that <I>2^^ is only about 3 times less intense than 
^^^^ in field h2. Those measurements do, however, agree 

(2) 

with our prediction that is several times more in- 
tense than ^^^^ in h2. 

VI. INELASTIC NEUTRON-SCATTERING 
MEASUREMENTS 

In earlier work^^ with D' = 0^ we obtained conflict- 
ing estimates for the easy-axis anisotropy K based on 
the spectroscopic and neutron-scattering spectra. Be- 
cause the instrumental resolution is broader than AnS/a 
[9], inelastic neutron-scattering measurements at the AF 
Bragg point (27r/a)[0.5, 0.5, 0.5] average over a range of 
q that includes both cycloidal satellites at (27r/a)[0.5 ± 
(5,0.5,0.5 =F For D' = 0, the spectroscopic mode fre- 
quencies indicated that K ^ 0.002 but the inelastic- 
scattering spectra indicated that K ^ 0.004. 
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FIG. 6: (Color online) (a) The measured inelastic-scattering 
spectrum^' around 77 = and the predicted spectrum for (b) 
K = 0.0025, (c) 0.003, and (d) 0.0035 meV with So ranging 
from to 0.015. 



We now re-examine the spectrum x^\uj) for ^ 0. 
The upper left-hand corner of Fig. 6 plots the measured 
spectrumf^i^. The resolution-averaged intensity spec- 
trum is plotted versus uj in Figs.6(b-d) for three values 
of K and six values of Sq from to 0.015. The very 
low- frequency rise of X^\oj) due to <l>o at 7^ = (5 has been 
removed from both the measured and predicted spectra. 

Below 5 meV, the measured x^^oo) contains four peaks 

at 1.2, 2.4, 3.4, and 4.4 meV. The peaks at 1.2 and 2.4 

(1 2) 

meV are primarily caused by ' ^ and ^o- As shown 

(2) 

in Fig. 4 for = 0.015, the separation between ^-j^ and 

$^^^/^o increases as K exceeds 0.001 meV. Correspond- 
ingly, the gap in the predicted spectrum centered at 2 
meV widens with increasing K beyond 0.001 meV. 

(2) 

As shown in Fig.5(b), <l>i is slightly enhanced by Sq. 
But the resolution-averaged spectrum x"{^) ^^Iso involves 
nearby modes and shifts to lower frequencies with in- 
creasing 5'o. For 5'o = 0.015 and K = 0.0035 meV, the 
low- frequency peak lies at 1.2 meV. So based on this sin- 
gle peak, K ^ 0.0035 meV provides good agreement with 
both the spectroscopic and inelastic measurements. Al- 
though its intensity increases with So and it is more pro- 
nounced than in our previous work^^, the predicted low- 
frequency peak at 1.2 meV is still considerably weaker 
than the measured peak. 

For K = 0.0035 meV, the second peak lies at 2.5 
meV when 6*0 = but shifts down to 2.3 meV when 
^0 = 0.015. More problematically, the predicted spec- 
trum contains three peaks between 2 and 4 meV (al- 
though the third peak is suppressed with So) whereas the 
measured spectrum contains only two. For K = 0.0035 
meV and 5*0 = 0.015, there are no predicted SW excita- 
tions between 4 and 5 meV at 7^ = or Consequently, 
the observed peak at 4.4 meV is missing from our spec- 



trum, which falls off much more rapidly than the mea- 
sured x"{^) above 4 meV. Keep in mind, however, that 
the predicted shape of x"{^) sensitively depends on the 
resolution function used to perform the averaging. 



VII. CONCLUSION 

A primary motivation of this work was to see how well 
a microscopic model can describe the properties of one 
of the simplest and most technologically important mul- 
tiferroic materials. We have demonstrated that all four 
modes observed by THz and Raman spectroscopies in 
BiFeOa are predicted by a model that includes two DM 
interactions, one along responsible for the cycloid pe- 
riodicity and the other along responsible for its tilt 
of the cycloid out of the x'z' plane. Using reasonable 
values for the easy- axis anisotropy and the DM inter- 
actions, we obtain excellent agreement with the mea- 
sured mode frequencies. The parameters D = 0.11 meV, 
D' = 0.054 meV, and K = 0.0035 meV provide very 
good descriptions of both the spectroscopic and inelastic 
neutron-scattering measurements, thereby resolving an 
earlier disagreement^^. 

The spectroscopic modes evolve with the complexity 
of the cycloid. With a single DM interaction D = Dy\ 
the cycloid is coplanar and purely harmonic. For nonzero 
frequencies, the only spectroscopically- active mode is ^1 
((^|M^/|0) ^ 0, ((5|M^/|0) ^ 0), which coincides with 
the EM (((^IP^^'^IO) ^ 0). Easy-axis anisotropy K along 
z' distorts the coplanar cycloid and introduces higher 
even harmonics in the Hamiltonian H. The 2Q poten- 
/(^) f/^\]\/r .ln\ ^ n /Alpindi 



tial splits ^±1 into (((5|M^/|0) 7^ 0, ((5|P^y^|0) 7^ 0) 
and (((5|M^/|0) ^ 0); the 4Q potential splits ^±2 
into <l>2^^ and $2^\ Hybridized with ^0 by the 2Q po- 
tential, <l>2^^ {{S\Myf\0) ^ 0) becomes spectroscopically 
active. Finally, the DM interaction = D^z' tilts the 
non-coplanar cycloid out of the x'z' plane. Then, ^0 
(((5|M3,/|0) ^ 0) and <l>^^^ {{S\My^\0) ^ 0) are dynami- 
cally and spectroscopically activated by their hybridiza- 
tion with ' and <l>o, respectively. Thus, additional 
interactions modify the mode spectrum as more modes 
hybridize with <l>o and ' \ 

Several experiments indicate that the low- 
temperature, low-field cycloid of BiFeOs under- 
goes a transition at about 140 K or 10 T. In THz 
measurements^^, the low-frequency ^o/*^i^^ mode dis- 
appears above 120 K and the high-frequency ^^"^ mode 
disappears above 150 K. Nevertheless, the selection rules 
governing the ^^^'^^ modes do not change^^. In Raman 
measurements, all modes persist for all temperatures 
but their frequenciesi^ and intensities display kinks at 
about 140 K. Optical^ and electron-spin resonance^ 
measurements show anomalies at about 10 T with 
indications that the cycloidal phase above 10 T is the 
same as the one above 140 K. Recently, Nagel et ai [l5[ 
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found that the THz modes exhibit kinks at about 5.5 T. 
But the nature of these transitions and the difference 
between the two cycloidal phases remain unknown. 

With magnetic field along z', the Hamiltonian of 
Eq.(|2]) does not produce a transition between different cy- 
cloidal phases^. Therefore, the proposed model may be 
incomplete. Since D' is responsible for the low-frequency 
^o/*^i^^ mode, a sudden change in D' at 140 K or 10 
T would produce anomalies in its spectroscopic features. 
A jump in D' at 140 K would also produce a jump in 
the weak ferromagnetic moment Mq(T). We hope that 
future experimental and theoretical work will resolve this 
and other mysteries surrounding BiFeOs. 
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Appendix A: SW intensities 

This section describes how to evaluate the SW inten- 
sities and eigenvectors X, which are required in the next 
section to evaluate the spectroscopic matrix elements. 

The local reference frame for each spin on site i is 
defined in terms of the unitary matrix by = IP S^. 
For spin 



and Hi must vanish. 



While £^0 is the classical energy 



where Vn 



(2M)\ 



(A4) 



) is a 4M- 



dimensional vector and ^(q) is a 4M-dimensional matrix. 

(r) 

Boson operators aq with 1 < r < M = 222 reside on 
layer 1 of the unit cell while those with M + 1 < r < 2M 
reside on layer 2. The sublattice index r refers to sites on 
either layer with R-x' = [r]a/\/2 where [r] = mod(r, M). 

and ttq^^^ obey the commutation relations 



Smce ttq 

Xr) (s)h _ 



and [a^\a[^)] 



0, Vq and vj^ sat- 



isfy the commutation relation [vq, vL] = X(5q q/ where 



iV-( ^ 
— I -/ 



and J is the 2M-dimensional unit matrix. 
A diagonal form for H2 is given by 



(A5) 



H2 



(A6) 



where Wq = (aq ^ 



.(2M) ^(l)t 



(2M)U 



q — yucd , . . . , aq % alq' , . . . , alq ' ' ) aud the 

boson operators aq^'^ and aq^'^^ also obey canonical com- 
mutation relations. The 4M-dimensional matrix L^(q) 
is diagonal with real eigenvalues en(q) = a;n(q)/2 > 
(n 1,...,2M) and en(q) -cJn(q)/2 < (n = 
2M + 1, . . . , 4M). So for each q, there are 2M positive 
and 2M negative eigenvalues. The commutation rela- 
tions yield 



Ho 



^c..(q)|a(-)ta(-) + l|, 

n,k ^ 



(A7) 



S = S (sin cos (/), sin sin (/), cos 6) , 
the matrices L[ and [/~^ are given by 



U 



cos U COS ( 

— sin (j) 
sin 6 cos ( 



cos U sm ( 

cos (/) 
sin sin c 




(Al) 



(A2) 



which identifies a;n(q) as the SW frequency for mode n 
with wavevector q. 

Vectors Wq and Vq are related by Wq = X(q) • Vq or 
Vq — X~^(q) ■ Wq, where the 4M-dimensional matrix X 
is normalized by X • TV • = TV. For fixed q. 



^(A,(q)-%6n(q))x^(q)=0, 



(A8) 



u- 



cos cos ( 
cos sin (f 
— s'mO 



— smc 
cos (j) 




sin cos ( 
sin sin 
cosO 



(A3) 



so that SU~^ • z = S. 

A Holstein-Primakoff transformation is used to express 
the local spin operators in terms of the bosons and 
a| with Siz = S — alai, 5'^+ = V2Sai, and Si- = \/26'a|. 
The Hamiltonian is then expanded in powers of l/V^ as 



where £(q) = L(q) • N. The inverse X"^ = TV • X"^ • TV is 
required to evaluate ((5|P^"^|0) and ((5|M|0). 

The wavevector Q and harmonic coefficients of the cy- 
cloid are obtained by minimizing Eq using the "trial" 
spin state provided by Eqs.(j2]|4]). If the spin angles on 
site r of layer 1 are Or and then the angles on layers 
1 and 2 are related by 6>^+m = Or -\- tt and ^r+M = — (/>r- 
We assume that (j)r = t and (pr+M = —t are independent 
of site position r on layers 1 and 2. 
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The spin-spin correlation function is defined by 



ttN J ^ 
^5(u;-u;„(q))5i";(q), 



iq-(Ri-Rj) 



(A9) 



where the final expression assumes that the SWs are un- 
damped. The inelastic neutron-scattering cross section 



a, (3 

= ^(l - {qjqr)s{u; - c^„(q))5i"j(q), (AlO) 

which only involves the diagonal matrix elements of 
Saisi^.^^) (if there is a net moment, some off-diagonal 
matrix elements a ^ j3 are nonzero and antisymmetric). 
The diagonal SW intensities (q) are given by 



S 



2M 



(All) 



where 



+ +^t^«,'0^.+2M,n+2M(q)- (A12) 

Even in the absence of damping, the instrumental res- 
olution will broaden the delta functions in 6'(q, cj) in 
Eq. ([AlOp . The magnetic form factor for Fe^+ should also 
be included in 5'(q, cj). 



Appendix B: Spectroscopic matrix elements 

This section evaluates the matrix elements for the in- 
duced electric polarization P^"^ and the magnetic mo- 
ment M between the ground state |0) and an excited 
state \S) with a single magnon at the cycloidal wavevec- 
tor Q. 

Since P^^^ = 0, only the and components are 
considered. Expanded about equilibrium, P^?^ becomes 



pin 



{M 

-\-S[r-\-2]-\-M,y' — S\ 
M 



Sm Ur COS ( 



r-2]-\-M,y' 



-S[r-\-2],y' + S[r-2],y' 



E 



sm Ur sm ( 



=1 



^[r+2],x' ~ S[r-2],x' 



~^^[r-\-2]-\-M,x' ~ S[r-2]-\-M,x 



(Bl) 



After some work, we obtain the EM matrix element 
for SW mode n: 



M 



COs6>[^+2] Sm{(j)r - (/>[r+2]) + ^ COs(^r - ^[r+2]) 

\^ [r+2] ,n+2M ^ [r+2]+M,n+2M J ^ 

+ COs6>[^+2] Sm{(j)r - (/>[r+2]) - i COs{(j)r - 0[r+2]) 

f — \ 2iqoa 

\^ [r+2]+2M,n+2M ^ [r+2]+3M,n+2M J ^ 

- COsO[r-2] Sin(^^ - ^[r-2]) + ^COs(^r - ^[r-2]) 

f — X~^ \ -2iqoa 

\-^[r-2],n^2M ^[r-2]+M,n+2M )^ 

- COs6>[^_2] Sin(0^ - 0[r-2]) - iC0s{(j)r - (/>[r-2]) 



2]+2M,n+2M ^ [r-2]+3M,n+2M 



-2iqoa 



(B2) 



where qo = 27rS/a. 

Similarly, P^^"^ can be expanded as 



S[r-\-2],x' ~ S[r-2],x' 



{M 
Y COS Or 

~S[r-\-2]-\-M,x' + S[r-2]-\-M,x' 
M 

S[r-\-2],z' — ^[r-2], 



sin Or cos (j)r 



-S[r-\-2]-\-M,z' + S[r-2]-\-M,^ 



(B3) 



The EM matrix element z' for SW mode n is 

M 

iqoar 



{6\Pi^^\0) = XSJ-J2e^^ 

^^[^+2] + ^ COS Or sin (/)[r+2] 

\ 2iqoa 
^[r+2]+M,n+2M )^ 

7r,[r+2] - ^ COS Or siu ^[^+2] 



2+2M 



-xr] 



^2iqoa 



r+2]+3M,n+2M 
r,[r-2] + ^ COS Or siu ^[^+2] 



(^[r+2],rH 

\^^[r+2]+2M,n+2M 
(^[r-2],n^ 



2+2M 



2]+M,n+2M 
^r,[r+2] - ^ COS t/r Sm ^[^-2] 



Y-1 

^[r-2]+3M,n+2M 



(^[r-2]+2M,n+2M 

(B4) 
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where 

gr^s = COS Or COS Og COS (j)s + sin Or sin ^5 cos (j)r ' (B5) 

For K = 0.0035 meV and S'o = 0.015, <l>o has the smah 
matrix element ((5|P^?^|0) ~ 0.19, about 60 times smaher 

than ((^|P^^^|0) ^ 12.2 for 

The MR matrix element for SW mode n is much more 
simply given by 

2M 

{S\Ma\0) = V2SfiB ^e^^°^M 

sgn(M - r + 1/2) I^,^,^HQ), (B6) 



which uses 

giQ.R _ gigoaH gg^(^ _ ^ ^ ;l/2). (B7) 

Notice that Wr^ (q) also enters the SW intensity S^J (q) 
of Eq.jAlI]). While the SW intensity S^aJ{Q) is propor- 
tional to the sum of \Wr^ {Q)]"^ over r, the matrix ele- 
ment {S\Ma\0) is proportional to the Fourier transform 

ofW^^c}{Q) over r. 
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